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Abstract 

The jamming transition in the stochastic cellular automaton model (Nagel- 
Schreckenberg model) of highway traffic is analyzed in detail, by studying 
the relaxation time, a mapping to surface growth problems and the inves- 
tigation of correlation functions. Three different classes of behavior can be 
distinguished depending on the speed limit v max . For v max = 1 the model 
is closely related to KPZ class of surface growth. For 1 < v max < oo the 
relaxation time has a well defined peak at a density of cars p somewhat lower 
than position of the maximum in the fundamental diagram: This density can 
be identified with the jamming point. At the jamming point the properties 
of the correlations also change significantly. In the v ma x = oo limit the model 
undergoes a first order transition at p — > 0. It seems that in the relevant cases 
1 < v max < oo the jamming transition is under the influence of second order 
phase transition in the deterministic model and of the first order transition 

&t Vmax — OO. 
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I. INTRODUCTION 



The breakthrough in statistical physics due to the understanding of critical phenomena 
and the related development of methods has initiated vivid interdisciplinary research activity. 
The successful application of these methods demonstrates that the validity of the concepts 
to handle the problem of many interacting units reaches far beyond the traditional scope of 
statistical physics. 

Vehicular traffic represents, from the point of view of statistical physics, a far from 
equilibrium driven system where a combination of the Highway Code and individual driving 
strategies replace the usual physical interactions between particles. Of course, the approach 
of a physicist is quite different from that of a traffic engineer: We would like to model the 
typical behavior as simple as possible while the essence of some phenomena should remain 
unaltered though it is not our aim to describe specific traffic situations. Nevertheless, one 
can hope a twofold gain from these type of studies: First, they could help us to understand 
far from equilibrium systems and, second, we strongly believe that the concepts developed in 
physics can contribute to the understanding of the complex phenomena related to vehicular 
traffic. Obvious analogies to problems of physics like to kinetic theory or granular flow have 
motivated physicists to work in this field. In fact, natural scientists' contribution to the 
understanding of traffic problems has already a long history ]I|,0|. 

Numerous articles have been published in the last couple of years investigating discrete 
models of highway traffic flow f||-0. It was suggested, that very simple probabilistic models 
based on cellular automata can reproduce features of real traffic, including a supposed 
transition from low density laminar flow to a high density phase, where start-stop waves are 
dominant. The behavior of these simple models is very complex near this transition, and up 
to now, is still not well understood. 

The paper is organized as follows. In the next section we summarize the model and the 
method of simulation. In Section |nj the relaxation time analysis is presented. Section [TV 
deals with the relation of the traffic models to surface growth. The correlation functions are 
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analyzed in Section |V|. The special case of no speed limit is discussed in Section [VI|. The 
paper is finished with a conclusion. 



II. MODEL AND SIMULATION 

In this paper we consider the cellular automata introduced by Nagel and Schreckenberg 
to describe single lane traffic |§. The model consists of a one-dimensional array of L cells 
with periodic boundary conditions. Every cell has (v max + 2) states: it can be empty, or it 
can contain a car with velocity v — 0,1, ... , v max . The density of cars is p. We perform the 
following steps in parallel for all cars: 

— Acceleration: increase v by 1 if possible. (l' a ) 

— Deceleration: decrease v to avoid crash with the car in front. (1'b) 
— Randomization: decrease v by 1 with probability p if possible. (l' c ) 

— Movement: move forward v sites. (I'd) 

Despite its simplicity this model captures several aspects of highway traffic including the 
free flow - jamming transition. The parameter v max can be considered as a limit speed and a 
usual value for it is 5. We would like to stress the importance of the third step. The fact that 
the model uses braking noise characterized by the probability p is crucial. One could equally 
introduce random accelerations, but it can be shown, that these type of perturbations die 
out very quickly, see ||. An initial condition is needed to specify the model completely. 

A convenient way to investigate the model is to draw a diagram of flow versus density, the 
so-called fundamental diagram. It is a curve with a well defined maximum at a density p c . 
The occurrence of density waves is related to the non-linearity of the fundamental diagram 
and it is expected, that the jamming transition will occur somewhere near the maximum. 
At low densities the flow is "free" with very few waves due to fluctuations (l'c), which die 
out quickly, at high densities above the maximum start-stop waves dominate the system; it 
is in the "jammed" state. 

What is the nature of this transition? It has been suggested that, with increasing 
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density, the jammed regions grow in space and time and at the jamming point they form a 
1 + 1 dimensional interconnected infinite network. This would point to a percolation type 
picture for the process. However, the percolation transition depends on the geometry of 
the clusters, i.e., on the microscopic definition of jams which introduces some ambiguity. 
Moreover, it was demonstrated in [|TU[ that appropriate choice of the parameters leads, for 
any reasonable jam definition, to a percolation transition at densities much higher than the 
region where jams become dominant. In order to resolve this problem and to avoid the 
ambiguity related to the definition of the jams, in [10| another definition of the jamming 



transition point was given without explicitly referring to the density waves. It was found 
that relaxation time of the average velocity has a maximum at a well defined density p p 
and the value of the maximum increases with growing system sizes. This phenomenon was 



interpreted in [|10j as critical slowing down and the critical point in the infinite size limit was 
identified with the jamming transition which turned out to be somewhat below the density 
at the maximum in the fundamental diagram (p c > p p ). 

The interpretation of the jamming transition as a second order non-equilibrium phase 
transition rises the question of the order parameter. One suggestion has been [D]] to take 
the density of jammed cars as order parameter which shows a rapid increase numerically at 
the same density where the relaxation time has its maximum. Besides the above mentioned 
problem with the definition of the jam this "order parameter" has the disadvantage that, 
due to the random braking events, there are small jams for any nonzero car density, i.e., the 
order parameter would not entirely vanish even below the transition. 

The present paper is devoted to a detailed investigation of the jamming transition. We 
have studied the model with three different values of v max namely 1,2 and oo. Because of 
the boundary conditions the v max = oo case means v max = L. It is widely believed that the 
Vmax = 1 and the 1 < v max < oo cases differ qualitatively [6] and we have found that the 
case without speed limit shows interesting and unusual features. 

With the upper choice of the v max parameters we can compare the three classes of models. 

The stationary state of the v max = 1 model is analytically exactly solvable f6j. In 



addition the model has a car-hole symmetry which connects the stationary state at p with the 
stationary state at (1 — p). In the v max > 1 cases this symmetry does not exist because every 
v > 1 step corresponds to the common movement of several holes. Due to its relationship 
with other non-equilibrium models the v max = 1 model is not expected to exhibit a "phase 
transition" at a specific density while this is assumed for v max > 1. 

For this reason we made simulations for the different cases and methods. We used p = 0.5 
for the braking probability throughout the paper. For the v max = 1 and the v max = oo cases 
we used an algorithm with position and velocity coding, storing the velocities and positions 
of the cars in two arrays. For the 

Umax — 2 CclSG WG used a multi-spin coding algorithm where 
we stored the lane v max + 1 times having L bits for every velocity according to the system 
size. In the lane labeled by v the z-th bit is 1 when there is a car with velocity v at site % and 
bit zero otherwise. This algorithm is faster because we can use bitwise operations. When 
calculating velocity correlation functions, we used the multi-spin coding algorithm only for 
reaching the stationary state and then switched to the position coding technique what was 
more appropriate for the calculations afterward. 



III. RELAXATION TIME ANALYSIS 

First we measured the average car velocity v = J2i 

function of time. As an 

initial condition we used uniformly distributed cars with Vi — for all cars and averaged 
over several runs. The performed simulations are listed in TABLE | and TABLE ||. 
We used the same r definition as introduced in flOf , namely: 



(mm(v*(t),( Voo )) - (v(t)))dt (1) 
Where v*(t) is the velocity-time function of cars without interaction and is the average 



velocity at the stationary state. As in |T(J we found a peak in r as a function of the density. 
For the v max = 2 case r m the maximum of this peak is at p p ~ 0.15 which is smaller than 
p c = 0.25 [FIG. 0,FIG. p|. The reason of the appearing negative relaxation times we see 
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at FIG. H] is that in this density region the average velocity function v(t) has a maximum 
after a fast increase for short times and it relaxes from this maximum to the stationary state 
[FIG. |3[]. Every timestep during this relaxation gives a negative contribution to r. 
Assuming that 

T m {L)~L z (2) 

we can fit z « 0.286 ± 0.012. On the other hand cr(L) the half-width of the r peaks is 
roughly the same. Therefore, in contrast to our earlier study with v max = 5, we have here 
no implication to assume a 

a{L) ~ IT* (3) 

scaling form and we could not do any finite size scaling. 

In the v max = 1 case p c = 0.5 as well known and p p « 0.35 [FIG. |]. The curves seem 
to be more similar than in the v max = 2 case and they seem to have the same form for 
different system sizes . The value of the z exponent is z = 0.27 ± 0.02. This value seems to 
be universal for all the v max parameters taking notice of [|10j too. However, in the v max = 1 
case this behavior is not bound to p p but it can be obtained in a broad range of the density. 

Looking for a different definition of a relaxation time we tried to fit a 

v(t) = Voo (l -b*T x *exp (--)) (4) 

function on the measured v(t) functions. Here Voo is the average velocity in the stationary 
state and b, x and r are fitting parameters. We began the fit after the first increasing regime 
which consists of 3-10 steps because then the car movements are nearly independent and 
therefore the v(t) graph is roughly a straight line. 

For every fit the value of the r exponent arose as a greater value than the corresponding 
system size. The exponent was vague to fit and had a great error because its effect falls in 
the regime of the fluctuations around foo(p) the stationary value of v(t). In the v max = 2 
case we could fit only in a restricted density region because of the shape of v(t) described 
above. 
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In the v max = 1 case the received values of the x exponent were the same within errors. 
Their value is 0.62 ± 0.09. This value is near to 2/3. This value can be interpreted in terms 
of a mapping to the KPZ surface growth problem (see the next Section). 

In the v m ax = 2 case the value of the x exponent increases from 0.6 to 1.5 in the measured 
density region [FIG. |5|. We have an upper density limit for the fits because of the emerging 
maximum in v(t) [FIG. [|. 



IV. RELATION TO SURFACE GROWTH 



With an appropriate conversion this traffic model can be transformed into a surface 



growth model [E|. The transformation is as follows. We go along the lane from left to right 
and order a slant line to every lattice site whose length is ^/2x lattice spacing and bevels 
with an angle a from the lane. At every step we order a line with a slope of a = 45° to a 
car and a line with a slope of a = —45° to a hole. We always continue this emerging zig-zag 
line at the end of it and we make a step upward or downward according to whether we find 
a hole or a car. After every update of the traffic model we get a surface in this model. The 
update of the traffic model determines the update of the surface growth model [FIG. ||. 

If the density p ^ 0.5 the surface has a density dependent average steepness of (2p — 1)\/2 
because the difference of the two ends of the surface is V2(2p — 1)L. 

In the v max = 1 case the surface growth picture corresponds to a deposition model in 
which we drop squares with their corner downward in every local valley with probability 
(1 — p) as we go from left to right. The squares have an edge of v^x lattice spacing. The 
deposition of a square corresponds to a step with velocity v — 1 in the traffic model. 

This simple deposition picture does not hold for the v max = 2 case because there are 
v = 2 steps also which would correspond to the deposition of a rectangle or correlated 
deposition of squares. 

In our simulations we measured the time dependence of the average surface width 

w 2 (t) = m,t)-h(t)Y) l (5) 



where h(i, t) is the height at site i, (.)$ means averaging over the sites and h(t) is the average 
height function 

h= (h(i,t))i. 

In a broad class of growth models w(t) behaves according to the Kardar-Parisi- Zhang theory 
13fl that is 



w ~ 



t> 3 



if t < U 



(6) 



and 



w ~ L c 



if 



t > U 



(?) 



The simple deposition model corresponding to the v v 



1 case is known to belong to this 



universality class [12 



As initial condition we used equally distributed cars. For the uniformly distributed 
initial condition we received large fluctuations and have not seen any regularity. In that 
case we saw the same behavior as in the system with equally distributed initial condition 
after relaxation. We subtracted the average steepness from the surface height in order to 
calculate w 2 (t). The performed simulations are listed in TABLE [TTI| and TABLE |TV| . 

In the v max = 1 case in the whole density region we see KPZ-like behavior [FIG. [7| and 
we received 2/3 = 0.63 ± 0.04 and 2a = 0.99 ± 0.06 which should be compared to the KPZ 
values (3 = 1/3 and a = 1/2. 

This KPZ-like behavior explains the 2/3-close value of the x exponent received from the 
relaxation time analysis (§) because the velocity of the cars: 



v(t + l) = W ((h(i, t + 1) - h(i, t)) = ±(h(t + 1) - h{t)) 

Z 8=1 A 



(8) 



that is, 



-v{t) 



dh 
~dt 



v s (L,t) . 



(9) 



According to []T4] : 



Av s (L, t) = v s (L, t) - v° s ~ t~ a± for t < , (10) 

where t>!? is v s (L,t) for £, L — > oo . The exponent corresponds to our x exponent and is 
in one dimension x = a± = 2/3. 

In the v max = 2 case the shape of the w 2 (t) curves is also KPZ-like. We measured 2a to 
be 2a = 1 ± 0.05. In contrast to the v max = 1 case the (3 exponent depends on the density 
as visible on [FIG. || and it is not uniform as in the v max = 1 case. 

V. CORRELATION FUNCTIONS 

Finally we studied the velocity correlation function in the space of car series with the 
definition 

C(i, t) = ((v(j, t') - v(t'))(v(j +i,t' + t)- v(t' + t))) j (11) 

or 

C'(i,t)^(v(j,t > )v(j+i,t > + t)} j (12) 

where v(j,t) is the velocity of the j-th car at time t and means averaging over all cars 
and t' is a time when the system is in the steady state. Being in the stationary state C(i, t) 
is independent of t'. 

If we consider the function C(i, t) at a constant time t as a function of i, we find a peak 
centered at the value 

i(t) = max C(j,t) . (13) 
5 

The maximum of the peak decreases with increasing t [FIG. [|. i(0) = because that is the 
autocorrelation function but for greater times it appears at other car indexes. Considering 
the i(t) versus t graph we can fit a straight line on the i(t) values which lines steepness 
depends on the density. This steepness defines a velocity in the car series space with which 
the peak spreads. 
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Except for the peak the C(i,t) function fluctuates around a constant value and it forms 
a plateau. 

From the definition of C(i, t) it follows that the jammed cars do not contribute to C(i, t) 
since they are staying. The cars already accelerated near v max give contribution to the 
plateau of C(i,t). Those cars give contribution to the peak which accelerate in a correlated 
manner because a car accelerating makes place to the car behind to accelerate. These are 
the cars coming out of a jam. 

Therefore the longterm existence of the peak and its velocity characterizes the traffic 
jams in the system. 

In the v max = 1 case the Vi(p) velocity of the peak increases without any sign of criticality 
from zero to its value at p = 1 where v(p — * 1) = (1 — p). This value follows from the motion 
of one single hole [FIG. |1C|| . 

However, in the v max = 2 case we can see a bending at the V^(p) function at the density 
Pk = 0.125 where the function shows a steep increase [FIG. |TT|. This implies emerging traffic 
jams. The small V^p) values of the function below pk come from the small fugitive jams 
being in the system. It is tempting to interpret the bending in V^p) as the appearance 
of a new phase at pk, since the behavior of the peak velocity is reminiscent to an order 
parameter of second order phase transitions. However, the velocity does not go to zero at 
the bending point and - in contrast to what is expected for second order transitions - we 
could not observe any finite size scaling. 

VI. THE Vmax = oo CASE 

Our simulations have indicated the following: 

i) The behavior v max > 1 is different from v max = 1 

ii) For v max > 1 there is a jamming transition but a careful study of the different character- 
istics indicated - in contrast to what was suggested in |10[ - that there is no clear second 
order phase transition related. 
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iii) From the present and earlier studies 0,|T^,|T^] it is clear that the smaller p the more 
pronounced is the transition. 

iv) Comparing our results for v max = 2 with those with v max = 5 we realize that increasing 
v-max also sharpens the transition. 

The origin of iii) is clear: In the p = model there is a phase transition related to a 
singularity in the fundamental diagram. In order to see the origin of iv) we considered the 
case without speed limit, more precisely the v max = L case. 

The fundamental diagram of the v max = oo case shows different features than the cases 
discussed so far [FIG. |l2j. The main characteristics of its shape is a plateau value depending 
on p [FIG. |13j. The function decreases fast to this plateau. The height of the plateau is 
independent of the system size only its length changes with. 

Looking at the flow of the cars we can observe that at the densities where the current 
of the stationary state joo(p) is higher than the j p plateau value the cars tend to be equally 
distributed what is typical for the deterministic case ||. However, at the density regime 
of the plateau we see one jam which characterizes this density region [FIG. [TJ|]. For higher 
densities where joo(p) decreases with the density we see more than one jam. 

Calculating the ( |T2] ) correlation function we can notice that the foo(p) peak velocity in 
the plateau region corresponds to j p which also means that we have only one jam in the 
system for this densities. 

Assuming that the p t value where the plateau sets in depends on L like 

Pt ~ L- f (14) 

we got t ~ 0.5. Because this point represents the Nf = p t L number of cars being in the flow 
between the right and left side of the jam we can conclude that the density contribution of 
cars in the flow part scales as 

Pf ~ L ~ l ( 15 ) 

and therefore the density contribution of the cars being in the jam will be 
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Pi = p- Pf 



(16) 



Pj — > p when L — > oo. (17) 

Thus values greater than the j p plateau value of joo(p) (except the value at p — 0) can be 
considered as fluctuation effects. When p — > we have finite number of cars with infinite 
velocity (one car with velocity L — 1). In the thermodynamic limit 

joo(p -> 0) = 1 as L -> oo. (18) 

Therefore we can say that in the thermodynamic limit we have a phase transition at 
p = between the phases of joo(0) = 1 and no jams and of the j p value with one jam. 
Increasing p the plateau ceases and there are more than one jams in the system. 

This result contradicts the result of the mean field theory described in [[J where the 
fundamental diagram starts in the p = j = point with infinite slope. 

For smaller braking probabilities than p = 0.001 the plateau ceases and only a break in 
the graph remains. 



VII. CONCLUSIONS 

According to our comparison of the three different models we can conclude to the fol- 
lowing. As we summarized at the beginning of Section 6 there is a difference between the 
behavior of the models v max = 1 and v max > 1. In contrast with the v max = 1 model in 
the v max = 2 case we see a kind of transition but it is not a strict second order transition 
it seems to be a crossover transition. We found that the density at which the transition 
takes place is smaller then the p c = 0.19 value mentioned in |TB[. On the grounds of the 



obtained results for v max = oo where we found a first order transition we can conclude the 
the increase of the v max parameter makes the transition more striking. Our results indicate 
that the behavior of the models labeled by different p and v max parameters is guided by the 
transition points of the models with parameters p = or v max = oo. 
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TABLES 



TABLE I. Performed simulations for the relaxation time analysis for v max = 1 

System size Number of runs 

L = 2000 10000 

L = 4000 5000 

L = 8000 2000 



TABLE II. Performed simulations for the relaxation time analysis for v max = 2 
System size Number of runs 

L = 2048 10000 
L = 4096 5000 
L = 16384 1000 
L = 32768 500 



TABLE III. Performed simulations for the measuring of w 2 (t) for v max = 1 
System size Number of runs 

L = 1024 500 
L = 2048 200 
L = 4096 50 



TABLE IV. Performed simulations for the measuring of w 2 (t) for v max = 2 
System size Number of runs 

L = 512 20000 
L = 1024 10000 
L = 2048 5000 
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FIGURES 

FIG. 1. Relaxation times versus density measured for the v max = 2 case for different system 
sizes. The maximum of the peaks is at the density p p ~ 0.15 and they have nearly the same 
halfwidths. 

FIG. 2. Fundamental diagram for the v max = 2 model measured for L = 2048. The maximum 
of the diagram is at p c = 0.25. 

FIG. 3. Average velocity versus time graph for v max — 2, L = 2048, p = 0.3125. 

FIG. 4. Relaxation times versus density graph measured for v max = 1- We see similar curves 
for different system sizes. 

FIG. 5. The x exponent versus density graph for v max = 2. 

FIG. 6. Determination of the corresponding surface growth model to the traffic model. 

FIG. 7. The square of the average surface widths versus time graph for the v max = 1 case 
for different densities. The system size is L = 2048 and the curves correspond to the following 
densities, p = 0.125 (N = 256), p = 0.25 (N = 512), p = 0.5 (N = 1024), p = 0.75 (TV = 1536), 
p = 0.96875 (N = 1984). On the log-log plot the the KPZ-like behavior is apparent. 

FIG. 8. 2(3 exponent measured in the v max = 2 model for different densities. 

FIG. 9. Correlation times versus car index for different times. The concerned system is 
L = 256, p = 0.5, v max = 2. 

FIG. 10. The peak velocity Vi(p) as a function of the density for v max — 1. 

FIG. 11. The peak velocity V 2 (p) as a function of the density for v max — 2. 
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FIG. 12. Fundamental diagram for v max = oo for L = 1000. 

FIG. 13. The plateau values for different p braking probabilities in the v max = oo model. 

FIG. 14. The one-jam phase of the v max = oo model. The horizontal direction means the 
lattice sites the vertical direction the consecutive timesteps. 
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